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We consider the thermal expansion, change of sound velocity with pressure and temperature, and the Poisson 
ratio of lattices which have rigid units (polyhedra very large stiffness to change in bond-length and to bond- 
angle variations) connected to other such units through relatively compressible polyhedra. We show that in such 
a model, the potential energy for rotations of the rigid units can occur only as a function of the combination ©i = 
{Oi — (V X Ui)/2), where 6i are the orthogonal rotation angles of the rigid unit i and Ui is its displacement. 
Given such new invariants in the theory of elasticity and the hierarchy of force constants of the model, a negative 
thermal expansion coefficient as well as a decrease in the elastic constants of the solid with temperature and 
pressure is shown to follow. These are consistent with the observations. 

PACS numbers: 65.40.De, 62.20.de, 62.20.dj 
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I. INTRODUCTION 

Most solids expand on raising temperature. There has been 
a large interest recently in materials that do the opposite. A 
particular example is ZrW208 in which the volume decreases 
linearly with temperature from about 10 K to about 1000 
K. The thermal expansion coefficient a = {\/V)dV/dT sa 
—2.6 X 1Q~^K~^ and is isotropic over this range'"^. Other 
examples are discussed by Tao and Sleight"*. These solids ap- 
pear to share the feature that they contain 'rigid units' . By this 
is meant that the complicated lattice structure of such solids 
is composed of corner or edge sharing polyhedra which are 
very stiff towards internal bond-stretching or bond-bending 
compared to the stiffness for change of the relative angles of 
the polyhedra and to other elastic forces which determine the 
sound-velocities. 

It was proposed very early^'^ that the dynamics of rigid 
units must be of primary importance in negative thermal ex- 
pansion in such solids. One of the first models proposed is 
one in which each polyhedra is connected to others by three or 
more polyhedra. In this case, rotation of any polyhedra must 
be accompanied by the rotation of the whole solid. Rotation of 
individual polyhedra, which was the idea behind the negative 
thermal expansion is then not a thermodynamic variable. Per- 
haps, to circumvent this problem, a model was introduced^ in 
which the atoms connecting the polyhedra was split in to two 
with a spring in between them. The analysis of such a model 
was carried out by standard methods of evaluating the dynam- 
ical matrix with some assumed potential functions. This gave 
what were termed "rigid unit modes" throughout the Brillouin 
zone which appear not to agree with the spectra measured by 
neutron scattering^ or with the temperature dependence of the 
specific heat in ZrW20g'* which appear anomalous only at low 
temperatures. 

A variant of the model^, represented in two dimensions by 
Fig. 1 of Ref. |9] for such crystals was proposed. The es- 
sential feature of this model is that the polyhedra are of two 
kinds, those that are three or more-fold connected and those 
that are only two fold connected. In such a model the two fold 
connected squares can be replaced by rigid rods. The general- 
ization of this to three dimensions is also given in Ref. [9] and 




Figure 1 : The unit cell in our problem. In the ground state at T = 
the distance between the neighboring squares is a (lattice constant). 
The displacement of the center of each square i is u^, and its rotation 
from the equilibrium is expressed through angle 9i. The rods 1q con- 
nect neighboring squares and they are very rigid. The other rigidity 
is due to the varying angle between rods and squares. 



has the same properties as the two-dimensional model. The 
polyhedra in ZrW208, for example satisfy such properties. 
The two dimensional model is the same as that represented 
by Fig. [T] here, with the springs between squares replaced by 
fixed rods. 

In Fig.l of Ref. |9] the squares do not change their shape 
of size and the only degree of freedom in the unit-cell is the 
angle 9. The area (volume in 3d) occupied by the solid is 
the sum of the area (volume) occupied by the solid squares 
(polyhedra) and the area (volume) between them. The latter 
are bounded by fixed perimeters (areas). Under this condition, 
the area (volume) is maximum when the perimeter (surface) 
is most symmetric. Any distortion from the condition of equi- 
librium leads to change in 9 with a change in area (volume) 
cx {59^. Since the thermal average {{59^) cx T for T > wq, 
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the bond-bending energy, a decrease in volume with tempera- 
ture proportional to temperature is to be expected. 

The model of Fig. 1 of Ref. i^, though adequate to describe 
the physics of thermal expansion, is however not a good repre- 
sentation of a solid. The only degree of freedom is the rotation 
angle 9. The non-Unearity of the problem leads to an interest- 
ing mathematical problem but there there are no elastic modes 
in such models which all solids must have. A simple modifi- 
cation of the model removes this serious shortcoming and in 
fact makes the problem solvable much more easily. 

The modification is to replace the rigid sticks or the rigid 
two-fold connected polyhedra with sticks or polyhedra with 
finite compressive stiffness to obtain the model represented in 
Fig.[T] It is important that this stiffness be large compared to 
the energy to change 0. As we will show, this hierarchy of 
stiffnesses — very large (so as to be taken to be infinity or 
rigidly constrained) for all internal co-ordinates of the four- 
fold connected squares compared to that of the spring which 
is in turn is large compared to the stiffness for changing 6 pre- 
serves the feature which gives negative a in the earlier model. 
It has the virtue of having finite elastic constants required of 
a solid so that sound velocity as a function of volume or pres- 
sure or temperature can also be calculated and compared with 
experiments. Recent experiments show that solids with a < 
also have a decrease in elastic constants both with temperature 
and with pressure''"'. We will also be able to explain this be- 
havior We will also show that these properties are consistent 
with the general relation between the Griineisen parameter 7 
which gives the change in vibration frequencies with volume 
and a. 

The two dimensional model is exhibited in Fig.[T] Just as in 
the previous paper^, a three dimensional model with the same 
properties can be shown to have similar properties. Instead of 
just one rotational degree, they have three rotational degrees 
of freedom per unit-cell. But the primary feature, the new in- 
variant in the energy, exhibited in Eq. ^ and its coupling to 
other degrees of freedom, has similar properties in the there- 
dimensional model. Our intent in this paper is only to exhibit 
the special features introduced by the hierarchy of the stiff- 
nesses in solids with rigid units and not to calculate the elastic 
properties quantitiatively. 

We should add that the physical features of the model are 
similar to those in the model of previous investigators^"^. Our 
formulation of the problem reveals new general physical and 
mathematical features of the problem. Some numerical work 
has also been done on a model similar to that treated here^. 



Ui = (Ri — Ri°'') = {ufx + u\y) is the displacement of the 
center of the mass of the square labelled iij) from its equi- 
librium position at T = 0. The two extra degrees of freedom 
compared to the earlier model ensure elastic properties of the 
model, the Oi degree of freedom provides an optic mode which 
is crucial as before for the anomalous thermal properties for 
T comparable or larger than its energy. 

Let us focus on a unit cell. The geometric position of the 
squares and springs are shown in the schematic figure Fig. [1] 
For convenience, we consider a single square at position i, and 
find its free energy. The summation over the entire lattice is 
then assumed. 

The two neighboring corners of squares i and i + x are 
spaced apart by vector 



1,: 



in 



(1) 



Unit vectors f,; = icosOi^smOi) and = (— sin0i,cos0i) 
point along diagonals of square i, and ^ is the ratio of square 
diagonal to the square (i.e., lattice) spacing a. The vector l^.a 
corresponds to the rod between squares i and i + x. Analo- 
gously, for the other three corners/rods we have 



= -ax + Uj_s - Uj + a|(fi + fi_£), (2) 



= —ay + u. 



- Uj + a|(si + Si^y). 



The nominal length of each of the rods at T = is = 
a(l-0- 

In this model there are two contributions to the potential 
energy. There is a rotational stiffness due to bending between 
the rods and a square; this is parametrized by an angle /^i ^ 
defined through (a = ±x or ±y) 



cos /ii,±j 

COS ^li^±y 



llli.itxil 



\h.±y\ 



(3) 



The other stiffness is due to the tendency of the rods to pre- 
serve their nominal length l'^^^ and it is proportional to the 
elastic constant of the rods g that is taken to be large com- 
pared to the rotational stiffness. The total potential energy is 



V = 



^K{1 - cos fii 



f(l|l..o| 



(4) 



II. EQUATIONS OF MOTION AND NORMAL MODES 

Let us count the degrees of freedom per unit-cell in the 
model. There are 8 total degrees of freedom per unit-cell 
and five rigid constraints which fix the lengths and the inter- 
nal angles of the squares. That leaves 3 dynamical degrees 
of freedom per unit-cell. (In the earlier model^, the length 
of the rods was also fixed introducing two more constraints 
thus leaving only one degree of freedom, 9 per unit-cell.) The 
three degrees of freedom may be taken to be {ui,9i} where 



The potential energy Eq. (|4|l needs to be augmented with the 
kinetic term. Assuming each square has mass AI, the kinetic 
energy reads 



(5) 



The second term, the rotational energy, is proportional to the 
moment of inertia M{a^)'^J. The constant J is a dimen- 
sionless number depending on the mass distribution in each 
square. For example, the uniform mass distribution yields 
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J = 1/3, while mass concentrated at the square corners re- 
sults in J = 1. 

In Appendix A, we calculate the potential energy for long 
wave-length vibrations up to the fourth order in the displace- 
ments u and rotations 0.. The harmonic part of the potential 
energy, Eq.©, at long wavelengths is 



r\2 ,1/2 I 2 \ I 2 



(6) 



where n = 2K/{a^g{l — ^)^) is a dimensionless parameter. 
Also, 



e = (6*- X u 



(7) 



The second and the third terms in Eq. (|6]) are the familiar terms 
from the theory of elasticity, where we also introduced strain 
fieldsii 



dxi dx-i 



(8) 



The first term in Eq. ^ is unique to the model. A term in 
energy proportional to ( V x u)'^ is forbidden because a global 
rotation of the solid should cost no energy. However, a term 
proportional to 8^, as derived in the appendix, is obviously 
allowed. In fact the interaction expanded to any order in the 
displacement can depend on 9 and V x u only in the com- 
binations 8, as explicitly shown in the appendix to cubic and 
quartic order. This is the basic new feature of the vibrations of 
solids with rigid units. All the interesting new results of this 
paper can be traced to this feature. 

Arranging the degrees of freedom into a column ~ 
{6, , u^) the equations of motion to harmonic order can be 
concisely written as 



Mdiag[Ja2^^ 1, 1] + Myip = 0. 



In the long wavelength limit, we have 



Mv 



It is easy to verify that 




KlQy 









(9) 



(10) 



-iqy/2 




(11) 



Then the above vector is an eigenvector and we can rewrite 

t 

^:^e + iqyU^/2-iq^uy/2^Q (12) 




The corresponding eigenvalue is lo^ = 2a^ g k / {M J) = 
jj^^a ■ There are also two gapless elastic modes (linear com- 
binations of y) required of a two dimensional lattice with 
eigenvalues ' 



a g 

23792;, In- 



setting qx = q cos ( 
7^ can be written as 



and qy = q sin 4>, the dispersion up to 



(^3 r = 



AM 



[2 + K± 



(13) 



V(2 - 2k + k2) + 2(1 - K.) COS' 



2gK 



1 + ii^^±i±2Ma2^2 



(14) 



AK 



1 + ?(-W2JO a2 2 



In order to obtain the q^ term in Eq. ( fT4b we need to use the 
subleading long-wavelength expansion of the potential energy 
given in Appendix B. Up to this order, the dispersion of the op- 
tical mode, which is the mode of the 8 degree of freedom Eq. 
( fT4l ) is isotropic and does not depend on the spring constant g. 

The elastic modes have a velocity anisotropy consistent 
with the underlying lattice symmetry, and correspond to a 
square lattice with leading order elastic constants 



5, 



0, cl 



44 



5K 

2 



K 



(15) 



It should be noted that we have constructed a very simple 
model to highlight the new features due to the dynamics of 
rigid units and kept only the bare minimum of rigidities rele- 
vant for the anomalous thermal properties. Because we chose 
only nearest neighbor central force constants = 0. We 
can therefore make comparison to experiments only for the 
elastic constants Cn and C44. 

We will consider corrections to these elastic constants at 
finite temperature and with changing pressure in Sec. IV. 

The Poisson ratio of the model follows from the effective 
elastic constants Eq. ( fTSl l. For C°2 = the Poisson ratio is^ 



1^20 



- 2C°4) sin' 2(j) 
C°^{1 + cos2 2(1)) + 2CO4 sin' 2(f> 

(1 - k) sin' 2(j) 
(1 + cos2 2(j)) + Ksin' 241' 



(16) 



When K < 1, as it is in our case, the Poisson ratio is strictly 
positive. Conversely, when k > 1 it is always negative. So, 
it appears that not all solids with constraints have auxetic be- 
havior (i.e., < 0)'"'''"*. 

We can compare this spectrum of this model with the rigid 
constraint model. If we assume the angle 9 is small and the 
difference of displacement is high order than angles (u^+q — 
Ui) = 0{9f) for a — ii, ±y, we can approximately solve 



the constraint 1 1 L 



a(l - 0- We find 



V V 



(1-0 
U-0 



{2~Oi0f+9f^,) + 2^9.9.±i 



(17) 



i2-0{ef+9f^^) + 2^9,9,±y 



(18) 
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These solutions correspond to the solution Set II from Ref. 1 
Plugging these into the expression for the angles, we have 



1 — cos /ii_Q 



1 



4(1-^)2 

+2{-2 + ^)^e,e^+a 



{[2 + {-2 + OS.K0^ 



f)2 



J 

(19) 



Then we find the equation of motion for 9 is 



AK 



The resulting dispersion is 



K 



1 



(21) 



We compare this result to the optical mode Eq. ( fT4l i. The gaps 
are the same, however, the long wavelength dispersion differs 
slightly. We notice, however, that if we were to take first the 
limit g ^ oo, and k — !■ of the harmonic terms in a controlled 
way such that their product K = Kga^{l — /2 remains 
finite, the dispersions coincide. 



III. CHANGE IN VOLUME WITH TEMPERATURE 

So far we have only considered the harmonic approxima- 
tion for the new model. At this level, the thermal expansion 
must be zero. The change of volume 6V is proportional to the 
divergence of the displacement field. In harmonic approxima- 
tion: 



SV/V = (V • u) 



= 0. 



(22) 



Here and after the summation over indices is assumed in term 

In order to get non-zero 5V, we must consider anharmonic 
terms in our Hamiltonian which we find from the expansion 
ofEq. © given in the Appendix. The cubic anharmonic term 



(3) 



pot 



dr 



Index a my be either x or y. Eq. ( l24b is equivalent to the Eqs. 
(fTTj-fTsTl derived earlier from consideration of the constraints 
alone. 

Integrating Eq. ( |24] | and taking its average we arrive to 



(20) 



9i 



2(1-0 



9a (e^). 



(24) 



(v-u)(r) = - 



(25) 



The thermal average of 8^ is given below, and the classical 
limit is 



Afa2^V[(iw„)2-(c^3(q))2] 



1 



dq 



2Ma2^2jw3(q) 



coth 



/3w3(q) 



2a'^gK 



(26) 



In the last step we assumed that the massive optical mode is 
only weakly dispersive. The thermal expansion coefficient ob- 
tained by combining Eqs. (|25] | and (|26] | is 



2aV(l-e) 



fcB^(l-0 
AK 



(27) 



2(1 -^) {^'^^"^'^^ " 2^9^^ - Uyy){u^y) + 

^{Uxyf{u^^) " (1 " [{Uxxf + [Uyyf] }. (23) 

If we were to keep all the anharmonic terms the ther- 
mal expansion would depend on the thermal averages {Q'^)t, 
{{uij)Q)T, and {{uij){uki))T- We know, however, that in the 
physical limit that we are interested in (k, <C 1), the first aver- 
age is of order 0(k)^^, while the others are of order 0(k)". 
Thus they are smaller by a factor of k, which allows us to 
concentrate on the first term in Eq. (l23T l. For the same rea- 
son, terms proportional to n were omitted in Eq. ( |23] | and will 
be omitted hereafter in all other calculations. Adding this an- 
harmonic term to the Lagrangian the equations of motion (at 
a; = 0) become 



The physics of the negative thermal expansion coefficient 
is the same as in the earlier model and discussed in the In- 
troduction. The usual anharmonicity of the longitudinal and 
the transverse acoustic modes cannot change this because in 
the hierarchy of the assumed stiffnesses, such modes have 
low occupation compared to the rotation modes of the rigid 
units because their energy is higher than W3, Eq. ( fT4] i over 
most of the Brillouin zone. The thermal contraction (as well 
as the change in sound velocity with temperature derived be- 
low) comes from the first anharmonic term in Eq. ( |23] ) which 
in turn can be traced to the fundamental constraint, Eq. ([T]| 
which links the rotation freedom of the rigid units which costs 
low energy to the motion of the degrees of freedom which 
determine the volume of the lattice. 



IV. CHANGE IN VELOCITY OF SOUND WITH 
TEMPERATURE, PRESSURE OR VOLUME 

The change in the elastic constants with temperature at 
fixed pressure of ZrW20s^ is similar to those of the normal 
solids which expand on increasing temperature, i.e., they de- 
crease on increasing temperature. At a fixed temperature the 
elastic constants also decrease with pressure. It is worth dis- 
cussing in general that the two properties are mutually consis- 
tent and consistent with the relation of the Griineisen parame- 
ter 7 and the thermal expansion coefficient. 

The relation 



duj_ _ {doj/dP)T 
W^^ ~ {dV/dP)T 



(28) 
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gives that since {dV/dP)T < for stability, decrease in the 
vibration frequencies w with pressure must be accompanied 
by an increase in uj with volume. = (— f^)T is therefore 
negative. The thermal expansion coefficient is shown in gen- 
eral to be proportional to 7 defined to be an average over all 
cj of 7tj , on thermodynamic grounds^^. It thus follows that 
the the observed behavior of the change in volume with tem- 
perature, and the change in elastic constants with pressure are 
mutually consistent. However, we proceed to show this ex- 
plicitly by calculating both the change in the elastic constants 
with temperature and with pressure. 

Consider the cubic anharmonic terms in Eq. (|23] |. As seen 
above, a change in temperature leads to a nonzero expectation 
value for strains at finite temperature 



(29) 



Therefore term such as Uxx, and Uyy in Eq. ( |23] | can be re- 
placed by their expectation value. This leads to corrections to 
the harmonic part of the potentiaU^. This correction is linear 
inT: 



(30) 



3^ [K.)2 + K,)2]}(V.u)(T). 

Small correction due to k are ignored in Eq. ( [30l l. The first 
term of Eq. ( |30l l yields small corrections (proportional to T) 
to the mass of the optical mode and can therefore be neglected 
in our work. Only the second and the third term contribute to 
renormalization of the elastic constants Eq. ( fTST l. 

In addition to Eq. dSOl l. there is another term of the same 
order in temperature affecting the dispersion of elastic modes. 
It is derrived from quartic anharmonic term and the fact that 
the average {0'^)t also depends on temperature. Expanding 
potential Eq. (|4]i up to the forth order and keeping only 0(k)° 
terms, we find 



(4) 



pot 



/ dr 



4(1 -ey 



}.(31) 



Here we have introduced a parameter a to fix the relative mag- 
nitude of the quartic to the cubic anharmonicity. In the simple 
model that we have introduced, this parameter is necessary to 
get the right variation of all the elastic constants with tempera- 
ture as well as pressure. Terms of kind Q{du)^ and (du)'^ lead 
to higher order corrections and have been dropped. Substitut- 
ing the thermal average (6x)t in Eq- <EB we obtain another 
correction to the harmonic potential 



4(1-0^ 
-(l + 0[(".x)' + Ky)']}- (32) 



Same as in Eq. ( [30] l, the first term only affects the optical mode 
mass and is therefore neglected. 



Combining the two corrections Eqs. dSOl l and ( |32] | to poten- 
tial, the temperature dependence of effective elastic constants 
is found 

l^jp ^ ' 2(1 -e) ^''^ 



fda 



44 



V dT 



(l-3a)C 



(34) 



where a is the (negative) thermal expansion coefficient found 
in Eq. ( l27T i. 

The pressure dependence of the elastic constants at a con- 
stant temperature is derived in an analogous fashion. 



dCA 



dP 



(3^a)-(3 + a)g 

' 2(1^0 ^^^^ 

(1-3q)C .... 

-9 ^ _ ^ ^T, (36) 



where kt is the isothermal compressibility measuring change 
in volume proportional to applied pressure; {5V)/V = 
—kt{6P). Within the present model 



(V-u)(P) 



2 2 
-P = — P. 



C'n+C',, 9 



(37) 



The change of the elastic constants with pressure and tem- 
perature are consistent with the thermodynamic requirement 
discussed in Eq. 



V. COMPARISON WITH NORMAL SOLIDS WITH a > 

It is worthwhile discussing how the differences from solids 
with normal thermal expansion coefficient a > come about. 
We may consider the following (schematic) long wavelength 
Hamiltonian for them. 



H = KE+ {l/2)duKdu + gaSuaucfu. 



(38) 



KE is the kinetic energy, k is the matrix of harmonic force 
constants and ga the cubic-anharmonic tensor. An expectation 
value (dudu) = vo^ (T) is finite at all temperatures so that we 
may decouple the last term and minimize to get that there is a 
change in volume which may be schematically represented by 



{du) = -gsw^T)/^. 



(39) 



Solids usually have 93 < 0. This ensures that a in- 
crease(decrease) of interatomic distance decreases (increases) 
interatomic interaction energy relative to the harmonic case. 
This then leads to a positive thermal expansion. This is to be 
contrasted with the g derived above in the model with rigid 
units from the expansion of potential energy subject to the 
constraint Eq. In this case 53 > 0, as also required by the 
general argument given in the Introduction. 

Given a finite (du), we may decouple the second term in 
Eq. ( [38] l also in the alternative way to find an effective stiffness 
renormalization: 



Keff = Ko + 2g3{du) ^ K-2g^w {T)/k 



(40) 
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We see that the correction to the effective stiffness does not de- 
pend on the sign of (73; we always have a decrease in the elas- 
tic constants with temperature due to the cubic anharmonic 
term since the mean square displacement (T) must always 
increase with temperature consistent with entropy always in- 
creasing with temperature. However the quartic term acts in 
an opposite direction in general and to get elastic softening 
their effect should be smaller than that of the cubic anhar- 
monicity. As we see below to get this feature of the experi- 
ments we have to introduce the parameter a ^ 1. 

VI. COMPARISON WITH THE EXPERIMENTS 

Zirconium tungstate is about the most studied materials 
with negative expansion coefficient. We now check whether 
our calculations of the thermal expansion coeffcient Eq. dZTl i. 
and the temperature (Eqs. ( [33] l and ( [34l i) and pressure (Eqs. 
( |35] l and ( [36] l) corrections to the elastic constants are in quali- 
tative accord with the experiments. 

In our model a negative thermal expansion coefficient a, 
linear in temperature occurs for temperatures above the energy 
of the rotation optic mode. In experiments, a oc —T above 
about 10 K. We can estimate the parameter K « lOOKelvin 
using an effective mass of 0(500) atomic units, to get w 
lOK. 

The thermal expansion coefficient a = —2.6 x lO^^K^^. 
Our result yields this order of magnitude with K « 10^ K for 
^ w 0.9. 

The relative decrease of the elastic con- 
stant {l/C^i)dCn{T)/dT is 0(10-3) and of 
(1/044)9044 (T)/9r is an order of magnitude smaller—. 
Similarly the decrease of O44 with pressure is an order 
of magnitude smaller than that of On—. To get the ratio 
between these quantities, we need to use a « 1/4 in Eqs. (|34] | 
anddlSll. 

With these values, we calculate {l/C^^)dCii{T)/dT of 



the right sign but about 5 times larger in magnitude. This 
indicates the perils of trying to get the experimental numbers 
through a model which is too simple in terms of details de- 
spite our introduction of the phenomenological parameter a. 
Part of the difficulty stems from our ignoring normal anhar- 
monic terms to concentrate on the new physics due to the rigid 
unit rotations. The usual anharmonic terms are always used as 
adjustable parameters anyway. 



VII. CONCLUSIONS 

We have constructed a model of a solid with negative ther- 
mal expansion coefficient. The solid consists of rigid units 
interconnected by elastic units where the elastic stiffness is 
much higher than the next relevant stiffness related to the 
bond bending. The model has sensible elastic properties even 
though the optical mode is identical to that found earlier in 
fully constrained models, which did not have that virtue. A 
novel feature is the uncovering of a new invariant in the theory 
of elasticity in this class of models. This may be of relevance 
also in other contexts. 

We do obtain the correct qualitative features for the ther- 
mal expansion coefficient and the ratio of the change in the 
temperature and pressure dependence of the two elastic con- 
stants we have calculated. The quantitative agreement cannot 
be aspired for, although We obtain the thermal expansion co- 
efficient, which depends on the details less than the elastic 
constants, quite well quantitatively with reasonable parame- 
ters. We have had to introduce an ad-hoc parameter a, the ratio 
of the cubic to the quartic anharmonic parameter, to estimate 
the change of elastic constants with temperature and pressure 
quantitatively. Despite that this simple model gives the quan- 
titative value only to within a factor of about 5 although the 
relative change in the two elastic constants calculated is given 
much better. 



Appendix A: Expansion of the Potential Energy 

In this appendix section, we show the details of the calculation of the derivation of the expansion of the potential energy into 
harmonic and leading anharmonic terms. Since we assume the stiffness of the spring is much larger than the rotation, we will 
focus on the expansion of 



;(0))2 



(Al) 



In order to show the result in a rotational invariant form we introduce Q and strains: 



e = 61 - -V X u, 
2 



1 r 



(A2) 



The nonlinear term in the strains is important for proving that at each expansion order the potential energy is rotational 



invariant. It is straightforward to expand Eq. (lAli) to find the quadratic terms V2 = ^ {d^u^Y + {dyU^Y 



The cubic and 
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fourth order potential terms are given as follows. 



Vi = 



2(1-0 

2 

a g 
4(1^ 



^2qa ^ 2ee\~d^uy + dyu'-') + ^6''\{i + 2e)[(a.w^)' + [dyu'-'f] - 2e[(a,u")' + (dyuyf] 



Here • • • in V4 are 8(9u)^ and ((9u)^ type of terms which lead to higher order corrections. Plug 9 = Q + \ {dxU^ — dyU^) into 
V3, then rearrange terms, we find 



1^3 



2(1-0 



(A3) 



Now we want to write this in terms of strains. But we have to remember there are cubic order contributions from V2. we find 

V2 

In Eq. (IA4l) . 0^{du) denote the 4th or higher order terms of du. The second term of Eq. (IA4b should be combined with V3 



^ [{u,,)^ + (uyy)^] - ^ [id,u-)' + {dxu-){dxuyf + {dyuyf + ( 9, ) ( 9, ) ^ ] + 0\du) (A4) 



2(1-0 



^Q'^Uii - 2^Q{Uxx - Uyy)Uxy + Ui: 



2 l"'xx ^ "'yy 



+ o'^{du,e) 



(A5) 



Here 0'^{du, 8) are 4th or higher order terms in du and 8 which we ignored in V3 and ua = Uxx + Uyy 

We can go on to rewrite V4 in terms of 8 and strains. We should include the 4th order contribution from V2 and V3. But for 
V4, we only need terms proportional to 8^, 8"^ and 8^, thus the only correction terms we need is the following 



a'g^ 



■8^ 



{dxu^f + {dxuyy + {dyu^'' + (dyuy) 



x\2 



,V\2 



4(1-0 

which is coming from ^Q{dxU^ + dyuy) term in Eq. ( IA3b . Then after some algebraic calculation, we find 

a'g^ 



V4 



4(1-0^ 



^8'^ + 8' 



- (1 + K ) 



+ 0^(9w,8) 



(A6) 



These contributions are rewritten in Sec. 111. However, to compare with experiments, it is found necessary in Sec. Ill to introduce 
a phenomenological parameter a, which is the ratio of the quartic to the cubic anharmonic potentials. 



Appendix B: Subleading long wavelength expansion of the potential energy 



In this section we present the gradient expansion of potential energy Eq. (|4|i where the next order terms in lattice spacing a 
are not truncated. These terms are relevant for finding the weak dispersion of the optical mode in Eq. ( fT4l l and the comparison 
of our model to the fully constrained model at the end of Sec. II. For all other purposes, these terms can be truncated and the 
harmonic potential in Eq. ^ are sufficient. 

Using earlier defined phase space vector ip, the potential energy is written as 



Vpot — 



a^g 



2n 



V 



1 - 



dx 



6 



dx 



a'dt 



a'dt 



Uy 



Z '-14 

a d„ 
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6 



dl 



I 



In the constrained model, the limiting procedure where k 
remains finite yields the dispersion of Eq. (ISTT i. 
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